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Abstract 

Time-frequency analysis, such as the Gabor transform, plays an important role in many signal processing 
applications. The redundancy of such representations is often directly related to the computational load of any 
algorithm operating in the transform domain. To reduce complexity, it may be desirable to increase the time and 
frequency sampling intervals beyond the point where the transform is invertible, at the cost of an inevitable recovery 
error. In this paper we initiate the study of recovery procedures for non-invertible Gabor representations. We propose 
using fixed analysis and synthesis windows, chosen e.g., according to implementation constraints, and to process the 
Gabor coefficients prior to synthesis in order to shape the reconstructed signal. We develop three methods to tackle 
this problem. The first follows from the consistency requirement, namely that the recovered signal has the same Gabor 
representation as the input signal. The second, is based on the minimization of a worst-case error criterion. Last, we 
develop a recovery technique based on the assumption that the input signal lies in some subspace of Li- We show that 
for each of the criteria, the manipulation of the transform coefficients amounts to a 2D twisted convolution operation, 
which we show how to perform using a filter-bank. When the under-sampling factor is an integer, the processing 
reduces to standard 2D convolution. We provide simulation results to demonstrate the advantages and weaknesses of 
each of the algorithms. 

I. Introduction 

Time-frequency analysis has become a popular tool in signal processing. During the past three decades, it has been 
successfully used for noise suppression UJ, (2), blind source separation [3|, echo cancelation J4), 0, JSJ, relative 
transfer function identification [7 |, beamforming in reverberant environments (8), system identification in general 
AD, iflOl . and more. In algorithms based on time-frequency transforms such as the Gabor representation, there 
is often a tradeoff between performance and computational complexity, which can be controlled by adjusting the 
redundancy of the transform. The latter is determined by the product of the sampling intervals in time and frequency, 
which we denote by a and b respectively. Specifically, as a and b are increased, there are less coefficients per time 
unit for any given frequency range, and therefore the amount of computation needed to process the them decreases. 
This effect is especially notable in adaptive algorithms, where a directly affects the convergence rate. 

This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice, after which this version 
may no longer be accessible. 

This work was supported in part by the Israel Science Foundation under Grant no. 1081/07 and by the European Commission in the framework 
of the FP7 Network of Excellence in Wireless COMmunications NEWCOM++ (contract no. 216715). 

The authors are with the department of electrical engineering, Technion-Israel Institute of Technology, Haifa, Israel (phone: +972-4-8294682, 
fax: +972-4-8295757, e-mail: ewa.matusiak@univie.ac.at, tomermic@tx.technion.ac.il, yonina@ee.technion.ac.il). The third author is also with 
the University of Vienna. 



July 19, 2009 



DRAFT 



2 



The signal processing literature on Gabor-domain algorithms heavily relies on the fundamental requirement that 
any signal can be recovered from its coefficients in the transform domain. This requirement leads to the upper bound 
ab < 1. However, since the performance-complexity tradeoff is of continuous nature, it seems very restrictive to 
limit the discussion to this regime. Specifically, by increasing the sampling intervals beyond this bound we may 
further reduce the computational load of any algorithm operating in the Gabor domain. This benefit is obtained 
at the expense of additional deterioration in performance. It is important to note that when ab > 1, an additional 
source of performance degradation comes into play, which is the inherent reconstruction error. This is because in 
this regime, we can only guarantee perfect reconstruction for signals lying in certain subspaces of L 2 , as we show 
in this paper, but not for the entire space. Nevertheless, the resulting complexity reduction may be of greater value 
in some applications. 

In this paper, we explore reconstruction techniques for non-invertible Gabor transforms, namely in which ab > 1. 
The fact that in these cases perfect recovery cannot be guaranteed for every signal introduces extra flexibility 
in choosing the analysis and synthesis windows of the transform. Specifically, we address the case where both 
the analysis and synthesis windows of the transform are specified in advance. They can be chosen according to 
implementation considerations, for example finite support windows, or certain multiple-pole windows [ 1 1 1 admitting 
an efficient recursive implementation. Our goal, then, is to process the transform coefficients before reconstruction 
such that the recovered signal possesses certain desired properties. 

To tackle this problem, we borrow several approaches from the field of sampling theory, which has reached a 
high degree of maturity in recent years 11121 . lfl3ll . We begin by employing the consistency criterion in which the 
recovered signal f(t) is constructed such that its Gabor transform coincides with that of the original signal f(t) 
lfl4l . We then proceed to analyze a minimax strategy, where the reconstruction error || / — / || is minimized for the 
worst-case input f(t) 1131 . Both these approaches are prior- free in the sense that they do not make use of any 
special properties, or prior knowledge, that might be available on the signal. 

A prevalent prior in the sampling literature, is that the signal to be recovered lies in a shift invariant (SI) subspace 
of L 2 (see e.g., |[T6l . ifTTl . ifTHl . |[T9l . |20|, 1211 and references therein). In fact, signals and images encountered in 
many applications can be quite accurately modeled as belonging to some SI space lfl2l . lfl3ll . such as the space of 
bandlimited functions, the space of polynomial splines and more. Their widespread use can also be attributed to the 
link that subspace priors have with stationary stochastic processes l22ll . Il23ll . E4l . ||251 . which have been shown to 
constitute realistic priors for the behavior of natural images (26). In this paper, we generalize the STprior used in 
the sampling community to a richer type of subspaces of L 2 , which we term shift-and-modulation (SMI) invariant 
spaces. The third class of inverse Gabor techniques we consider, then, makes use of the SMI prior. We show that 
such a prior can lead to perfect recovery in some cases, given that the synthesis window is chosen according to 
the prior. For a fixed synthesis window, which is not matched to the prior, we show how to achieve the minimal 
possible reconstruction error for signals in the prior-space. 

In each of the three techniques we develop, the processing of the Gabor coefficients amounts to a 2D twisted- 
convolution ll27l with a certain kernel, which depends on the analysis and synthesis windows. We show that the 
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twisted-convolution operation can be interpreted in terms of a filter-bank. Furthermore, in the case of integer 
under-sampling (i.e., when ab is an integer), the resulting process reduces to a standard 2D convolution in the 
time-frequency domain. In these cases, we discuss situations in which the 2D convolution kernel is a separable 
function of time and frequency. This allows a significant reduction in computation, namely by implementing the 
2D convolution as two successive ID filtering operations along the time and frequency directions. 

The paper is organized as follows. Section [TT] is devoted to notation that will be used throughout the paper. In 
Section Hn] we derive conditions on the analysis and synthesis windows such that they generate Riesz bases for their 
span, which guarantees that the non-invertible Gabor representation is stable. In Section [TV] we review sampling 
and reconstruction schemes in shift-invariant (SI) spaces in order later to be able to generalize them to the Gabor 
transform using SMI spaces. Sections [VI I VII and IVIII constitute the central part of the paper, where in the first two 
we develop prior-free recovery procedures for Gabor transforms in the integer and rational under-sampling regimes 
respectively, and in the last we discuss SMI-prior recoveries. We devote Section I Villi to describing how twisted 
convolution can be realized as a filter-bank and also how to obtain the inverse of a sequence with respect to twisted 
convolution. Finally, in Section [IX] we demonstrate the methods we develop for the case in which both the analysis 
and synthesis are performed with Gaussian windows. 

II. Notation and Definitions 

We will be working throughout the paper with the Hilbert space of complex square integrable functions, denoted 
by L 2 (R), with inner product 

/oo 
f(t)W)dt for all f,geL 2 (R), (1) 
-oo 

where g(t) denotes the complex conjugate of g(t). The norm, induced by this inner product, is given by 

ll/f = </,/>• (2) 
The Fourier transform of / E L 2 (M) is defined as 

J7M= / f(t)e- 2 " u "dt. (3) 



For convenience, we will sometimes write / for Tj. 

In order to ensure stable recovery we focus on subspaces of L 2 (R) which are generated by frames or Riesz 
bases. A collection of elements {sk}k&z is a frame for its closed linear span if there exist constants A > and 
B < oo such that 

^||/|| 2 <^|(/,s fc )| 2 <5||/|| 2 for all /espaS{ Sfc }, (4) 

where span denotes the closed linear span of a set of vectors. The vectors {sfcjfcez form a Riesz basis if there 
exist A > and B < oo such that for all sequences c E i 2 

A\\c\\% < \\j2cks k f < B\\c\\% , (5) 
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where ||c||| a = J2kez\ Ck \ 2 * s ^ e sc l uare d ^ 2 -norm of c k . A direct consequence of the lower inequality is that the 
basis functions are linearly independent, which means that every function in spanjsfc} is uniquely specified by its 
coefficients c k . 

The fundamental building blocks of the Gabor representation are the so-called Gabor systems. To define a Gabor 
system, let a > and b > be such that ab — q/p with p and q relatively prime, and let T ak and M b i, for k,l S Z, 
be the translation and modulation operators given by 

T ak f(t) = f{t -ok); (6) 
M w /(t) = e™ bH f(t). (7) 

For s G L 2 (R), the Gabor system Q(s,a,b) is a collection {M b iT ak s(t) ; (fc,Z) G Z 2 }. The composition 

M bl T ak f(t) = e 2 * Mt f(t - ak), (8) 

which is a unitary operator, is called a time-frequency shift operator. Many technical details in time-frequency 
analysis are linked to the commutation law of the translation and modulation operators, namely 

M bl T ak = e 2 ™ bkl T ak M bl . (9) 

When p = 1, the time-frequency shift operators commute, i.e. Mu T ak = T ak M b i, because e 27riabkl = 1 for all 
k,l G Z. One consequence of the commutation rule, which we will use in our exposition, is the relation 

(/, M W _ b „T afc _ am /) - e ^"(l-n)rn (M bn T am f, M H T ak f) . (10) 

When p = 1 this becomes (/, M bl _ bn T ak _ am f) = {M bn T am f, M bl T ak f). 

For s G L 2 (]R), the collection Q(s,a,b) is a Riesz basis for its closed linear span if there exist bounds A > 
and B < oo such that 



<B\\c\\% ce£ 2 , (11) 



A\\c\\ 2 2 < I ^2 Ck,iM b iT ak s 
and is a frame when 

A\\f\\ 2 < J2 \(f, M biT ak s)\ 2 < B\\f\\ 2 for all / e spa5{M M T ofc s} . (12) 

A necessary condition for Q(s, a, b) to constitute a frame for L 2 (R) is that ab < 1, 11281 . Moreover, if Q(s, a, b) is 
a frame, then it is a Riesz basis for L 2 (R) if and only if ab = 1 ll28l . In this paper we focus on the regime ab > 1, 
where Q(s,a,b) does not necessarily span i 2 (M). 

With a Gabor system Q(s,a,b) we associate a synthesis operator (or reconstruction operator) S : t 2 (I?) —> 
L 2 (R), defined as 

Sc = Ck,iM u T ak s{t) for every c 6 £ 2 (Z 2 ). (13) 

fc,iez 

The conjugate S* : L 2 (R) — * £ 2 (Z 2 ) of S is called the analysis operator (or sampling operator), and is given by 

S*f={(f,M bl T ak s}} for every / G L 2 (M). (14) 
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(a) Analysis filter bank (b) Synthesis filter bank 



Fig. 1: Filter-bank representation of the Gabor transform (a) and of the inverse Gabor transform (b). 



IIL Stable Gabor representations 

The Gabor representation of a signal f(t) comprises the set of coefficients {ck,i}k.iez obtained by inner products 
with the elements of some Gabor system Q(s,a,b) [28|: 

c k ,l = (f,M u T ak s). (15) 

This process can be represented as an analysis filter-bank, as shown in Fig. da). Consequently, s(t) is referred to 
as the analysis window of the transform. If Q(s, a, b) constitutes a frame or Riesz basis for L 2 (M.), then there exists 
a function v(t) 6 L 2 (R) such that any f(t) 6 L 2 (]R) can be reconstructed from the coefficients {ck,i}k,iez using 
the formula 

f{t) = c k,iMbiT ak v{t). (16) 

k,l£Z 

The Gabor system Q(v,a,b) is the dual frame (Riesz basis) to Q(s,a,b). Consequently, the synthesis window v(t) 
is referred to as the dual of s(t). The recovery process can be represented as a synthesis filter-bank, as shown in 
Fig. Hfc). 

Generally, there is more than one dual window v(t). It can be shown that any function v(t) satisfying (y, Mi/ a Tk/bs) = 
SkSi is a dual window. The canonical dual window is given by v = Q _1 s, where Q is the frame operator associated 
to s(t), which is defined by Qf = J^k zez (/> ^blT a k s ) MuT a ks{t). There are several ways of finding an inverse 
of Q, namely by employing the Janssen representation of Q, through the Zak transform method or iteratively using 
one of several available efficient algorithms ll28l . 
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In this paper, we are interested in Gabor systems that do not necessarily span L 2 (R) but rather only a (Gabor) 
subspace. A Gabor space is the set V of all signals that can be expressed in the form ( TTol l with some norm-bounded 
sequence Ck.i- Since perfect recovery cannot be guaranteed for every signal in L 2 (R) in these situations, we have 
the freedom of choosing the analysis and synthesis windows according to implementation constraints. However, in 
order for the analysis and synthesis processes to be stable, we would still like to assure that the systems Q(s, a, b) 
and Q(v, a, b) form frames or Riesz bases for their span. In this section, we give several equivalent characterizations 
of windows v(t) and sampling intervals a and b such that the Gabor system Q(v, a, b) forms a Riesz basis. 

For tractability, we assume throughout the paper that a and b are positive constants such that ab = q/p, where p 
and q are relatively prime. Moreover, we will consider only Gabor spaces whose generators come from the so-called 
Feichtinger algebra Sq, which is defined by 



So = { / G L 2 (R) | || /|| So := f \(f,MM)\dxdu < oo }, 



(17) 



where ip(t) is a Gussian window. An important property of Sq is that if v(t) and s(t) are elements from Sq 
then {(v, Mi,iT a ks)}k,i£i, is an £ 1 (Z 2 ) sequence. Examples of functions in Sq are the Gaussian and B-splines 
of any order. The Feichtinger algebra is an extremely useful space of "good" window functions in the sense of 
time-frequency localization. Rigorous descriptions of Sq can be found in [28 1 and references therein. 

The first characterization of Gabor Riesz bases we consider is stated directly in terms of their generator v(t). It 
is a simple corollary of a result on Gabor frames for L 2 (R), see 



Proposition III.l. Let v(t) s Sq and ab = q/p with p and q relatively prime. The collection Q{v, a, b) is a Riesz 
basis for its closed linear span if and only if there exist constants A > and B < oo such that 

AI P <V(lu,x) < BI p for almost all (uj,x) £ R 2 , (18) 

where I p is the p x p identity matrix and V(u>,x) is a p x p matrix-valued function with entries given by 



V r ^,x) = ^^v(x-ar-?^y(^-as- l -y- Maku , r, s = 0, . . . ,p - 1. 
G(v,a, b) is an orthonormal basis if A = B = 1. 



(19) 



Proof: By the Ron-Shen duality principle [29], Q(v,a, b) is a Riesz basis (orthonormal basis) for its closed 
linear span if and only if the system Q(v, l/b, 1/a) is a frame (Parseval frame) for i 2 (R). The latter is a frame for 
i 2 (R) if and only if there exist constants A > and B < oo such that the so-called frame operator S vv , defined 
as S vv f{t) = J2kdez(-f> Ml /a T k/b v ) Mi /a T k / b v(t) satisfies 

Tl < Sw < x J > (2°) 
b b 

where / is the identity operator on L 2 (R). This means that S vv is bounded and invertible on L 2 (R). It was shown 
in 11281 that, since l/{ab) — p/q, the operator S vv satisfies (f20b if and only if ( TT8l is satisfied, which completes the 
proof. ■ 
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Note that lo is a frequency variable associated with the discrete-time variable k, and similarly x is a time variable 
associated with the discrete frequency index I. Another valuable observation is that that V{u>, x) is a (1/a, 1/6)- 
periodic function. Furthermore, it has been shown in [28| that V ryS (uj, x) is continuous. Therefore, the lower bound 
in ( f]~8b can be replaced by the requirement that det(V(u>,x)) > for all (u>,x) S [0, 1/a) x [0, 1/6). 

The next characterization we consider is in terms of the twisted convolution operator. Specifically, the Riesz basis 
condition implies that Q(v,a,b) is a Riesz basis for its closed linear span if and only if there exist constants A > 
and B < oo such that 

A (c, c) < (r vv tj c, c) < B (c, c) for all c £ ^ 2 (Z 2 ), (21) 
where the 2D cross-correlation sequence r vv [k, I] is defined as 

r vv [k, 1} = (v, M u T ak v) . (22) 
The operation \ represents the twisted convolution defined by 

(d\\c)[k,l]= J2 d m ,nCk- m ,i- n e- 2 ™ b «- n)m . (23) 

When p = 1, twisted convolution becomes standard convolution, because the exponential term in d23l equals 1 for 
all m,n,l € Z. Therefore, u(£) generates a Riesz basis if and only if the twisted convolution (standard convolution 
when p = 1) operator with kernel r OT [fe, Z] is bounded and invertible. Invertibility of this operator translates to the 
invertibility of the sequence r vv [k, I] with respect to \ (* respectively). Proposition IIII. 1 1 states then, that this twisted 
convolution operator is bounded and invertible if and only if the matrix-valued function V(ui,x) is bounded and 
invertible for all u and x. Explicitly finding the inverse of a sequence with respect to twisted convolution is not a 
trivial task. We will address the problem in Section IVIIII 

Our last representation follows from restating Proposition IIII. II using a different, but equivalent, matrix-valued 
function that involves the cross-correlation sequence r vv [k, I] defined earlier. This new representation was first 
introduced in ll30l to characterize the invertibility of general Gabor frame operators. 

Proposition III.2. [30 J The matrix-valued function V(ui,x) of \1% coincides with the matrix-valued function 
Q{tjJ, x) whose entries are given by 

* r , s (W, X) = J2 r »v [s-r+pk, l] e -^iablr e -2ni(blx+akw) ^ (24) 

and therefore Q(v,a,b) is a Riesz basis for its closed linear span if and only if there exist constants A > and 
B < oo such that 

AI p < x) < BI p for almost all (uj,x) eM 2 . (25) 
In the integer under-sampling case p = 1, €>(w, x) of d24l i reduces to the scalar function 

*(«,*) - J2 r vv [kAe~ 2 ^ Ux+akul) = (Fr vv )(LO,x), (26) 
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where Tr vv is the 2D discrete-time Fourier transform (DTFT) of r vv [k,l}. Therefore, in this case condition (|25l > 
reduces to 

A < (Tr vv )(uj, x) < B for almost all (uj, x) E M 2 (27) 

for some A > and B < oo. 

The ^-characterization is of particular interest in our context as it can be used to investigate any twisted 
convolution operation with a sequence h G 1 1 (Z 2 ). Indeed, it was shown in [31 1 that such an operation is invertible 
if and only if the matrix-valued function 

a;) = J2 h s - r+pkJ e- 27Tlabrl e- 2 ™ ( - blx+ak ^ (28) 

fc,Z6Z 

is invertible for all u> and x. In fact, in some sense the function $(u>, x) is to twisted convolution what the DTFT is 
for convolution. Specifically, we show in Appendix [A] that for two sequences Cf-i and dk,i having ^-representations 
<& d (ui,x) and & c (uj,x) respectively, the matrix-valued function $>( ci,d \ui,x) associated to the twisted convolution 
c t| d, can be expressed as 

$ {ct,d) (u;,x) = & d (u>,x)& c {u,x). (29) 

We conclude this section with the observation that having a Riesz basis for a Gabor space V, it is possible to 
construct many others using equivalent generating functions. 

Proposition III.3. Let G(v, a, b) be a Riesz basis for its closed linear span V and ab = q/p with p and q relatively 
prime. Let 

w(t) = J2 hk,iM bl T ak v(t), (30) 

where {hk.i} is a sequence of weights. Then G(w, a, b) is an equivalent Riesz basis for V if and only if there exist 
constants A > and B < oo such that the {p x p) -matrix-valued function & h (u;,x) of H28i satisfies 

AI P < $ h (uj,x)& h (uj,x) H < BI P for almost all (u,x) G R 2 , (31) 

where & h (uj,x) H denotes the conjugate transpose of <& h (uj, x). 

Proof: See Appendix IE1 ■ 
In the case of integer under-sampling (i.e., when p = 1), & h (u,x) becomes a scalar function, which is simply 
the 2D DTFT of h^.i. In this setting, condition OTb becomes 

A < \<S> h (uj,x)\ 2 < B for almost all (w, x) E R 2 . (32) 

IV. Sampling and Reconstruction in Shift- Invariant Spaces 

To address the recovery of a function f(t) from its non-invertible Gabor transform, we will harness several 
strategies which were initially developed in the context of sampling theory. Specifically, the last two decades have 
witnessed a substantial amount of research devoted to the problem of recovering a signal f(t) from the equidistant 
point-wise samples of its filtered version, using a predefined reconstruction filter lfl2l . |fl~3], (32). As can be seen 
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/(*) 



s(-t) 



t = ak 



v(t) - -^/(t) 



£ S(t-ak) 

k— — oo 

(a) Sampling (b) Reconstruction 

Fig. 2: Sampling (a) and reconstruction (b) with given filters. 



in Fig. |2j the sampling stage in this setting, corresponds to the central branch in the analysis filter-bank of the 
Gabor transform shown in Fig. [TJ a )- Thus, the time-frequency plane is sampled in this scenario only on the lattice 
{{ak, 0)}/cgz. Similarly, the reconstruction process of Fig.|2]can be identified with the central branch of the synthesis 
filter-bank of Fig. [Tib). 

The main goal in this setting is to produce a set of expansion coefficients {d k } by processing the samples {ct}, 
such that the recovered signal f(t) possesses certain desired properties. In this section we briefly review three 
methods for tackling this problem, each based on a different design criterion. For more detailed explanations and 
a review of other methods, we refer the reader to ifPfl . Il33l . ifTSI . iTPl . Il32l . In the following sections, we will 
extend these results to the Gabor scenario. 

For simplicity, we assume here that a — 1. The reconstruction process of Fig. [2] can be written in operator 
notation as / = Vd, where V : £ 2 — » L 2 (R) is the synthesis operator associated with the functions {v(t — fc)}fcez, 
defined as 

Vd = d kv(t -k) = Y^ d k T k v(t) for every d G £ 2 (Z). (33) 

fcez fcez 

Similarly, since c k — (f(t), s(t — k)), the sequence of samples {c k } are obtained by applying the synthesis operator 
S*, which is the conjugate of the analysis operator S associated with the functions {s(t — k)}k^z- 

S*f = {(f(t),s(t-k))} = {(f,T k s)} for every / G L 2 (R) . (34) 

We will refer to S — span{u(t — k)} and V = span{u(t — k)} as the sampling and reconstruction spaces respectively. 
Spaces of this type are called shift-invariant (SI). 

As in the Gabor transform, we will focus on cases where the sets of functions {s(t — n)} and {v(t — n)} constitute 
Riesz bases for their span. Then, both the sampling and reconstruction are stable procedures. It is well known ||34| 
that the functions {v(t — n)} form a Riesz basis for their span V if and only if there exist constants A > and 
B < oo such that 

A < 4> vv (u) < B for almost all u G E, (35) 

where 

<f> vv (u) = -L Vl^-fc)! 2 (36) 

fcez 
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is the DTFT of the cross-correlation sequence 

r vv [n] = («(*), »(t-n)> = (v(t) * v{-tj) (n), (37) 

and v(lo) is the Fourier transform of v(t). In other words, {v(t — n)} is a Riesz basis if and only if the sequence 
r vv [n] is bounded and invertible in the convolution algebra 1 1 (Z, *). In particular, the functions {v(t — n)} form an 
orthonormal basis if and only if A = B = 1. Notice the analogy with condition d25l l (and d27l i in the case p = 1), 
which was developed for Gabor systems. 

A. Consistent reconstruction 

Perhaps the most intuitive demand from the recovered signal f(t) is that it would produce the same sequence of 
samples {ck} were it re-injected to the sampling device of figure [2{ a), namely 

(/(t), s(t - fc)) = c fc = </(t), s(t - k)) (38) 

for all k 6 Z. This consistency requirement was first introduced in T\A\ in the context of sampling in SI spaces 
and then generalized to arbitrary spaces in [35), [33 1. There, it was shown that consistent reconstruction is possible 
under the direct-sum condition <S © V = L 2 (R), where © denotes a sum of two subspaces that intersect only at 
the zero vector. This means that S and V are disjoint and together span the space L 2 (R). 
In the SI setting, the direct-sum condition translates into the simple requirement that 



\(f>sv(u)\ > A, for almost all w e M (39) 

for some positive constant A, where 

cp sv {Lo) = - k) (40) 

fcez 

is the DTFT of the cross-correlation sequence r s „[n] = (s(t),v(t — n)). Under this condition, reconstruction can 
be obtained by convolving the sample sequence {ck} with the filter /i con , whose DTFT is given by [14|, [36], [37] 

^con(^) = — -^r, (41) 

to obtain the sequence of expansion coefficients {dk}- 

If S and V are two arbitrary subspaces of L 2 (K) satisfying 5 © V = L 2 (R) (namely not necessarily SI spaces), 
spanned by the functions {s n (t)} and {v n (t)} respectively, then the sequence of expansion coefficients d can be 
obtained by applying the the operator 

H con = (S^)- 1 (42) 

on the sequence of samples c, where S and V are the synthesis operators associated with {s n (t)} and {v n (t)} 
respectively. The direct-sum requirement guarantees that S*V : i 2 — > £ 2 is continuously invertible. In the next 
sections, we will use this latter characterization to develop a consistent reconstruction procedure for non-invertible 
Gabor transforms. 
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B. Minimax regret reconstruction 

A drawback of the consistency approach is that the fact that f(t) and f(t) yield the same samples does not 
necessarily imply that f(t) is close to f(t). Indeed, for a signal f(t) not in V, the norm of the resulting reconstruction 
error f(t) — f(t) can be arbitrarily large, if S is nearly orthogonal to V. 

To directly control the reconstruction error, it is important to notice that f(t) is restricted to lie in V by 
construction. Therefore, the best possible recovery is the orthogonal projection of f(t) onto V, namely / = Pyf, 
a fact that follows from the projection theorem [38|. This solution cannot be generated in general, because we do 
not know f(t) but rather only the sequence of samples {c^} it produced. The difference between the squared-norm 
error of any recovery f(t) and the smallest possible error, which is ||/ — Pvf\\ 2 = ll-Py^/ll 2 ' i s called the regret 
||39l . The regret depends in general on /(t) and therefore generally cannot be minimized uniformly for all f(t). 
Instead, the authors in |[T5l proposed minimizing the worst-case regret over all bounded-norm signals f(t) that are 
consistent with the given samples, which results in the problem 

m inrnax||/-/|| 2 - \\P v ^f\\ 2 , (43) 
fev f eB 

where B — {/ : S* f — c, ||/|| < L} is the set of feasible signals. 

It was shown in [15] that the minimax-regret reconstruction can be obtained by filtering the samples with the 
filter h mx whose DTFT is given by 

#mxM = -—- r -, (44) 

where (f>yg(u), 0ss(cl>) an d ^v( w ) are as in (l40l > with the corresponding substitution of the generators v(t) and 
s(t). Note that the solution is independent of the bound L appearing in the definition of B. 

If the sampling and reconstruction functions form Riesz bases for arbitrary spaces S and V (not necessarily SI), 
then the sequence of expansion coefficients d can be obtained by applying the operator 

H m = (V*V)- 1 S*V(S*S)- 1 (45) 

on the sequence of samples c. The operators V*V and S*S are guaranteed to be continuously invertible due to 
the Riesz basis assumption. This more general characterization will be used in the next sections to develop a 
minimax-regret recovery method for non-invertible Gabor transforms. 

C. Subspace-prior reconstruction 

The consistent reconstruction approach leads to perfect recovery for input signals that lie in the reconstruction 
space V 04). The minimax-regret method, on the other hand, leads to the best possible approximation f = Pyf 
for signals f(t) lying in the sampling space S 03]. Therefore, the two methods can be thought of as emerging 
from the prior that f(t) lies in a certain subspace W of L 2 (R), where W = V in the consistent strategy and 
W = S in the minimax-regret approach. In practice, though, it is often desirable to choose the sampling and 
reconstruction spaces according to implementation constraints and not to reflect our prior knowledge on the typical 
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signals entering our sampling device. Thus, commonly neither constitutes a subspace prior W, which is good in 
the sense that ||/ — fyv/H is small for most signals in our application. 

A generalization of these two methods results by assuming that / € W where W = span{u>(i — k)} for a 
generator w(t), which may be different than s(t) and v(t). If the subspace W satisfies the direct-sum condition 
S 1 - © W = L 2 (JR), then the solution / = P v f can be generated by filtering the sample sequence Ck with lfT31 

#subM - -—^j t— r, (46) 

where (w), 0sw(w) and are as in d40b with the appropriate substitution of v(t), s(t), and w(t). 

For arbitrary sampling, reconstruction and prior subspaces S, V and W (i.e., not necessarily SI), the coefficient 
sequence d can be obtained by applying the transformation 

H suh = {V*V)- 1 V*W(S*W)- 1 (47) 

on the sample sequence c, where 14^ is the synthesis operator associated to the prior functions {w n (t)}. This general 
formulation will be used in the next sections to derive a subspace-prior recovery technique for non-invertible Gabor 
transforms. 

V. Integer under-sampling 

In this section we address the problem of recovering a signal f(t) from its non-invertible Gabor transform 
coefficients {ck,i}, given by (fl5l l. using a pre-specified synthesis window v(t). We focus on prior-free approaches 
that do not take into account any knowledge on the signal f(t). Specifically, here we employ the consistency and 
minimax-regret methods discussed in the previous section to the Gabor scenario. To emphasize the commonalities 
with respect to the SI sampling case, and to retain simplicity, we begin the discussion with the case of integer 
under-sampling (p = 1). In the next section we generalize the results to arbitrary p. 

A. Consistent synthesis 

In the Gabor transform, the sampling (analysis) space S is spanned by the Gabor system Q(s,a,b) and the 
reconstruction (synthesis) space V is the span of Q(v, a, b). As discussed in Section HV-AI consistent reconstruction 
is possible if S ®A = L 2 (M.). In the case of SI spaces, this direct-sum condition translates to the requirement that 
the cross-correlation sequence {(s(i), v(t — has an inverse in the convolution algebra £ 1 (Z 2 , *). A similar 

condition is true in the setting of Gabor spaces. 

The next proposition characterizes the class of pairs of analysis and synthesis windows satisfying the direct-sum 
requirement in the integer under-sampling regime. 

Proposition V.l. Assume that Q(s, a, b) and Q{v, a, b) are Riesz sequences that span the spaces S and V respectively, 
and ab = q G N. Then S © V = L 2 (R) if and only if the function $> sv (uj,x), defined as 

= Y, r sv [k,l]e^ Ux+ak ^ = (Fr sv )(u,x), (48) 
feJez 
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is nonzero for all (u),x) € [0, 1/a) x [0, 1/6). Here, 

r sv [k,l} = (v,M bn T am s) (49) 
is the Gabor transform of the synthesis window v(t). 

Proof: It was shown in [33), for general Hilbert spaces, that if S and V are spanned by Riesz bases Q(s, a, b) 
and Q(v, a, b) respectively, then S © V = L 2 (R) if and only if the operator S*V is continuously invertible on I 2 . 
Here, S* and V are the analysis and synthesis operators associated with Q(s,a,b) and Q(v,a,b), respectively. By 
definition, for any sequence c £ £ 2 (Z 2 ) 

(S*Vc)[k,l] = { 

\rn,n£Z / 

= ^ Cfe_ m ,/_ n (v, M& n T am s) 

= (r s „ * c)[M]. (50) 

Hence, the operator S"*^ is simply a 2D convolution operator with kernel r sv [k,l] = (v, Mb n T am s) and S*V is 
invertible if and only if r sv [k,l] is invertible in the convolution algebra ^ 1 (Z 2 ,*). As shown in Section [HI] this 
sequence has a representation defined by (l28l i. which is its 2D DTFT in the case p = 1. A sequence is 

invertible with respect to convolution if and only if its DTFT has no zeros. Therefore, r sv [k, I] is invertible if and 
only if <£> sv (oj, x) ^ implying that S 1 - © V = £ 2 (IR) if and only if $ st, (w, x) ^ 0. ■ 
Assuming that indeed S^ffi^l = L 2 (M), we know from Section llV-Al that to obtain a consistent recovery, we must 
apply the operator H con — (S*V)~ 1 on the coefficients {ck,i} prior to synthesis. In the proof of Proposition lV.il we 
showed that S*V is a 2D convolution operator with the kernel r sv [k,l] of j49l . Therefore, (S'*F)~ 1 corresponds 
to filtering the Gabor coefficients with the filter h con whose 2D DTFT is given by 

H con (uj,x) = — — -. (51) 

<£ sv (w, X) 

This filter is well defined by Proposition IV. 1 1 since we assumed that the spaces generated by s(t) and v(t) satisfy 
the direct sum condition. 

Observe that during the operations of analysis and pre-processing of the Gabor coefficients e^;, we in fact 
compute a dual Riesz basis for the reconstruction space V. In case the synthesis and analysis spaces are the same, 
namely S = V, we compute the orthogonal dual basis. However, when the spaces are different we compute a 
general (oblique) dual Riesz basis for V. 

Proposition V.2. Let G(s, a, b) and G(v, a, b) be Riesz sequences that span the spaces S and V respectively, where 
ab is an integer, and assume that S © V = L 2 (R). Then a dual Riesz basis for the space V is Q(g, a, b) with 

g(t) = J2 h C on[m,n]T- am M- bn s{t) e S. (52) 
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where h co „[m,n] is the inverse of r sv [k, I] with respect to *. 

Proof: Any signal in V can be recovered from the corrected coefficients dk,i — (/icon * c)[k, I] via f(t) — 
Sfe zez dk,iMbiT a kv(t), where c k> i is as in (T5[ . Therefore, we may view this sequence as the coefficients in a basis 
expansion. To obtain the corresponding basis we note that by combining the effects of the analysis window s(t) and 
the correction filter H con of dBTl l, the expansion coefficients can be equivalently expressed as d k .i = (/, M b i T ak g) 
where 

g(t) = ^ hcon[m, n]T- am M-bns(t) G S, (53) 



Indeed, 



(f,M bl T ak ,g) = / f,M bl T ak \ ^ h con {m,n]T^ am M- bn s 

\ \ m,n£Z 

= \ f> ^ h con {m,n}M b iT ak T- am M_ bn s\ 

\ m,n£Z / 

= (/, )^ h CO n[m, n]Mbl-bnT a k-amS ) 
\ m,nGZ / 

= h con [m, Tl] (/, Mu~bnTak-amS) 

^ ^ ^con[^i ^] C/c— ml— n 
m,n£2S 

= (hcon * c)[fc,Z] = d fe(i . (54) 

Therefore, any / G V can be written as 

/(*)= £ (f,M bl T ak g)M u T ak v(t). (55) 

It can be easily verified, by Proposition IIII. 31 that C/(<7, a, 6) is an equivalent Riesz basis for 5. Furthermore, it can 
be checked that 

(M u T ak v, M bn T am g) = S (56) 
implying that Q(g, a, b) is a dual Riesz basis to Q(v , a, 6). ■ 

B. Minimax regret synthesis 

We now wish to develop a minimax-regret recovery method, similar to the SI sampling case of Section IIV-BI 
Specifically, we would like to produce a recovery f(t) for which the worst-case regret ||/ — /|| 2 — ||Pyj./|| 2 
over all bounded-norm signals f(t) consistent with the given Gabor coefficients {c k ,i}, is minimal. As men- 
tioned in Section IIV-BI the minimax-regret reconstruction can be obtained by applying the operator H mx — 
(V*V)~ 1 S*V {S* Sy 1 on the Gabor coefficients c k .i prior to synthesis. 
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From Section IV- A| we know that when p = 1, the operators V*V, S*V and S*S correspond to 2D convolutions 
with the kernels r vv [k,l], r sv [k, I] and r ss [k, I] respectively, which are given by d49b with the appropriate substitution 
of s(t) and v(t). Therefore, the minimax-regret recovery is obtained by filtering the Gabor coefficients Ck.i with 
the 2D filter h mx , whose DTFT is given by 

& sv (u) x) 

H m x(u,x) = -— — -. (57) 

$ ss (w, X)9 VV (U), X) 

Here, $ sl, (w,x), $ ss (lu,x), and $™(w,x) are the 2D DTFTs of r sv [k,l], r ss [k,l] and r vv [k,l] respectively. This 
filter is well defined by Proposition IIII.2I since we assumed that s(t) and v(t) generate Riesz bases for their span. 

C. Efficient implementation 

As we have seen, the two reconstruction approaches discussed above are based on 2D filtering of the Gabor 
transform Ck,i prior to synthesis. A significant reduction in computation can be achieved in cases where the 2D 
correction filter is a separable function of k and I, namely when hk,i — UkVi for two sequences Uk and Vk- In 
these situations, one can first apply the ID filter u% on each of the rows of Cf-j (i.e., along the time direction), 
and then apply the ID filter V[ on each of the columns (along the frequency direction). If, for example, hf-i is a 
separable finite-impulse-response (FIR) filter with N x N nonzero coefficients, then direct application of it requires 
N 2 multiplications per output coefficient, whereas only 2N multiplications suffice when implementing it using two 
ID filtering operations. 

Separable correction filters emerge when the cross-correlation sequences involved are separable functions of k and 
I. One such example is the case where s(t) and v(t) are Gaussian windows with variances a 2 and a 2 respectively 
and aba 2 /(a 2 + a 2 ) is an integer (recall that we also require that ab be an integer). Then r sv [k, I], r ss [k, I], and 
r vv [k, I] are all separable functions of fc and I, so that both the consistent and the minimax-regret filters are separable. 
More details on non-invertible Gaussian-window Gabor transforms are given in Section [IX] 

VI. Rational under-sampling 

We now generalize the results of the previous section to the case where the product ab is not an integer, but 
rather some rational number q/p with p and q relatively prime. The main difficulty here is the fact that the time- 
frequency shift operators do not commute when p ^ 1. Therefore, instead of standard convolution we will be faced 
with a twisted convolution, which is a noncommutative operation. This makes the techniques from Fourier theory 
inapplicable in a straightforward manner. 

A. Consistent synthesis 

Obtaining a reconstruction f(t), which is consistent with the Gabor representation Cjy of f(t), is possible if 
S 1 - © V = L 2 (R). As we have seen in Proposition IV. 11 in the integer under-sampling case p = 1 the direct sum 
condition translates to the requirement that the cross-correlation sequence r sv [k,l] be invertible in the convolution 
algebra £ 1 (Z 2 , *). In the setting of rational under-sampling, we have the following. 



July 19, 2009 



DRAFT 



16 



Proposition VI.l. Assume that Q(s,a,b) and Q(y,a,b) are Riesz sequences that span the spaces S and V 
respectively, and ab — q/p with p and q relatively prime. Then S © V — L 2 (R) ;/ and only if the (p x p)- 
matrix-valued function <f> sv (u),x) with entries defined as 

*™ n (w,x)= J2 r S v[n-m + pk,l}e- 2mablm e~ 2m(blx+akuj) m, n = 0, . . . ,p - 1. (58) 

is invertible for all (u),x) £ [0, 1/a) x [0,1/6), which is equivalent to det(3?(w, a;)) 7^ for all (u>,x). 

Proof: The proof is similar to the proof of Proposition lV.il Since s(t) and v(t) generate Riesz bases for S and 
V respectively, the condition S 1 - ffiV = L 2 (K) is satisfied if and only if the operator S*V is continuously invertible 
on £ 2 , where S* and V are the analysis and synthesis operators associated to Q(s, a, b) and Q(v, a, b) respectively. 
By definition, for any sequence c € ^ 2 (Z 2 ), we have 

{S*Vc)[k,l] = / c m , n M bn T amVl M bl T akS \ 

= Cm,n(v,M bl _ bn T ak _ am s)e- 2 ^ bl - bn ^ m 

Ck-rn.l-n \V, M hn l arn S) e v ; 

= (r M t|c)[fc,Z]. 
Therefore, 5*1^ is a twisted convolution operator with kernel 

r sv [k,l] = (v,M bl T ak s), (59) 

and S*V is invertible if and only if r sv [k,l] is invertible in the twisted convolution algebra £ X (Z 2 , \\). As shown 
in Section Hill this sequence has a representation 3> sv (uj, x) defined by d28l and so is invertible if and only if this 
matrix is invertible. Therefore, S 1 - © V = L 2 (M.) if and only if $> sv (lu, x) is invertible for all w and x. ■ 

Note that for p = 1, the above proposition reduces to Proposition IV. II When p ^ 1, we conclude from 
Proposition IVI. 1 1 that the direct sum condition translates to the requirement that r sv [k, 1} be invertible in the twisted 
convolution algebra, which can be checked by analyzing its ^-representation. An alternative method for checking 
weather r sv [k,l] is invertible with respect to \, is presented in Section IVHII It involves only the sequence r sv [k,l] 
without introducing the continuous variables to and x, making it more attractive in some cases. 

As in Section IV- Al to obtain a consistent recovery f(t), we have to apply the operator H con — (S 1 *!/)" 1 to the 
Gabor coefficients c k ,i- However, as opposed to the case p = 1, where H con was a standard convolution operator, 
here it corresponds to a twisted convolution operation. This is due to the fact that time-frequency shift operators 
do not commute for p ^ 1. Specifically, in the proof of Proposition IVI. II it was shown that S*V corresponds 
to twisted convolution with r sv [k, I]. Therefore, corresponds to twisted convolution with the sequence 

r~J~[k, I], which is the inverse of r sv [k, 1} in the twisted convolution algebra ^ 1 (Z 2 , \\). This inverse exists, since we 
assumed that the spaces generated by s(t) and v(t) satisfy the direct-sum condition, and it will be shown in the 
next section how to construct it. 
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One can write the twisted convolution relation between the Gabor transform Cf-j and the expansion coefficients 
dk,i in terms of their ^-representations. Specifically, since d — (S*V)~ 1 c, we have Ck,i = (r sv \\d)[k, 1} and therefore 

* c (w, a;) = x)$ sv (uj, x), (60) 

where $> c (u>,x), & d (tu,x) and $> sv (lo,x) are the p x p-matrix- valued "^-representations of the sequences Cfc j, dk,i 
and r sv [k,l] respectively, defined in d24l . Therefore, to obtain the sequence dk.i from the Gabor coefficients Ck,i, 
we apply a twisted convolution filter, whose $ function is 

H con (Lj,x) = $ sv (w,x)- 1 . (61) 

The twisted convolution operation can be modeled as a filter bank which is specified by the convolutional inverse 
of r sv [k,l], as we show in Section [Villi 

During the operations of sampling and pre-processing of the samples Ck,i we in fact compute a dual Riesz basis 
for the synthesis space V. If the synthesis and analysis spaces are the same, namely S = V, we compute the 
orthogonal dual basis. However, when the spaces are different we compute a general (oblique) dual Riesz basis for 
V. 

Proposition VI.2. Let G(s, a, b) and G(v, a, b) be Riesz sequences that span the spaces S and V respectively, and 
ab = q/p with p and q relatively prime. Assume that V = L 2 (R). Then a dual Riesz basis for the space V 
is Q(g,a, b) with 

g(t) = ^ h C0 „[m,n]T- am M- bn s(t) G S , (62) 

where h co „[m,n] is the inverse of r sv [k,l] with respect to \. 

Proof: Any signal in V, that has been sampled with the Riesz sequence Q{s, a, b) resulting in the coefficients 
Ck,i given by ( fT5l ). can be recovered from the corrected samples dk.i = (/icon t| c)[k, I], where h mn [k, I] — r~ v }[k, I] 
is the inverse of r sv [k,l] with respect to t|, via f(t) = J2 m n£zdk,iMuT a kv(t). This sequence may be viewed as 
the coefficients in a basis expansion. To obtain the corresponding basis we note that by combining the effects of 
the analysis window s(t) and the correction twisted-convolution filter h con [k,l], the expansion coefficients can be 
equivalently expressed as dk.i — (/, MuT a f.g) where 

g(t) = ^ heoa[rn,n]T- am M-bns(t) G S. (63) 

Indeed, 

(/, M u T ak g) = I f, M u T ak \ ^2 h con {m,n]T- am M- bn s 

\ \ m,nEZ 

= ( /, ^ /icon [TO, n]M b iT a kT- am M-bnS ) 
\ 771, n^Z / 

= (f,J2 /iconK n]e 2mab ^- m ^M bl ^ bn T ak ^ amS \ , (64) 

\ m,n£Z / 



July 19, 2009 



DRAFT 



18 



and using the linearity of the inner product, we have 

(f,M bl T ak g)= h mn [m,n]e- 2naHk - m ^ (f, M bl ^ bn T ak _ am s) 

= (h coa \[c)[k,l] = d k> i. (65) 
Therefore, any / G V can be written as 

/(*)= J2 (f,M bl T ak g)M u T ak v(t). (66) 

It can be easily verified, by Proposition IIII. 31 that Q(g, a, b) is an equivalent Riesz basis for S. Now, for it to be 
a dual Riesz basis to Q(v, a, b) we need to check that 

(MuT ak v, M bn T am g) = 5 (67) 

Indeed, 

(M bl T ak v,M bn T am g) = e ^ iab ^ k (v, M b(n _ T a(m _ fc)P ) 

= e 2 ™ b(l - n)k hoon[x,y] (v,M b(n _ l _ y) T a(m _ k _ x) s)e- 2 ' iab ( m - k -^ 

^2Triab(l — n)k \ ^ t r v r_ j j i — 2iriab(m — k—x)y 

= e K } h, CQn [x,y\r sv [m ~ k ~ x,n - I ~ y\e K IU 

x.y^L 



= e 2 * lab{l - n)k {h m ^r sv )[m -k,n-l} = <5 m _ fe <5„_ ; , 



where we used the fact that M bn T am g{t) = J2x, y ez h con[x,y]e 27rtab( - m x)y M b{n _ v) T a{m _ x) s{t) and that h con is 
the inverse of r sv with respect to t|. ■ 

B. Minimax regret synthesis 

Next, we develop a minimax-regret reconstruction method for non-invertible Gabor transforms with rational 
under-sampling. Our goal here, as in Section IV-BI is to minimize the worst case regret max/ e g{||/ — /|| 2 — 
||Pyj_/|| 2 }, where B is the set of bounded-norm signals whose Gabor coefficients coincide with c k j. As dis- 
cussed in Section IIV-BI The recovery / attaining the minimum can be obtained by applying the operator H mx = 
{V*V)~ 1 S*V{S* on the Gabor coefficients c k: i prior to synthesis. However, as opposed to the integer 
under-sampling case discussed in Section IV-BI where V*V, S*V, and S*S were convolution operators, here 
they correspond to twisted convolutions with r vv [k,l], r sv [k,l] and r ss [k,l] respectively. Therefore, to obtain the 
sequence d k i, we apply a twisted convolution filter on the Gabor coefficients c k i, whose impulse response is 

h mx [k,l] = (r-^r sv ^)[k,l]- (68) 
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Here, r v ^ [k, 1} and r ss 1 [k, I] are the inverses of r vv [k, I) and r ss [k, 1} with respect to tj. Consequently, the $ function 
of the minimax-regret filter is given by 

H lm {u,x) = S^z)" 1 *^, (69) 

where $> ss (us, x), <fr sv (uj, x), and <& vv (ui, x) are the ^-representations of r ss [k,l], r sv [k, I] and r vv [k, I] respectively. 

C. Extension to symplectic lattices 

Throughout the current and previous sections, we considered a special type of sampling points in the time- 
frequency plane, called separable lattices A = aZ x WL. However, with the help of metaplectic operators, these 
results carry over to the more general class of lattices, called symplectic lattices. A lattice A s C R 2 is called 
symplectic, if one can write A s = T>A where A is a separable lattice and V G GL2(R), meaning it is an invertible 
2x2 matrix with determinant 1 |28|. To every T> G GL2(M.) there corresponds a unitary operator n(D), called 
metaplectic, acting on i 2 (M). One can show that a Gabor system on a symplectic lattice is unitarily equivalent to a 
Gabor system on a separable lattice under /j(£>), that is Q(g, A s ) is a frame/Riesz basis if and only if Q(fi(T>)~ 1 g, A) 
is a frame/Riesz basis, and 

g(g,A s ) = n(V)G([i(V)- 1 g,A). (70) 

Therefore, instead of considering a representation of f(t) in span{(7A}AeA s one can look at the representation of 
f{t) in span{/i(2?) _1 5A}AeA- For more details see |28|. 

VII. SUBSPACE-PRIOR SYNTHESIS 

In the previous two sections we attempted to recover a signal from its non-invertible Gabor representations 
without using any prior knowledge on the signal. When such knowledge is available, it can significantly reduce the 
reconstruction error and in some cases even lead to perfect recovery. A common prior in sampling theory is that 
the signal to be recovered lies in some SI subspace of L 2 , namely that it can be written as 

f(t) = J2 d k T ak w(t) = d k w(t - ak) (71) 

fcez fcez 

with some norm-bounded sequence {dk} and some window w(t). This model can quite accurately describe many 
types of natural signals, which exhibit a certain degree of smoothness. For example, the class of bandlimited signals 
is the SI space generated by the sine window. The class of splines of degree N also follows this description with 
w(t) being the B-spline function of degree N. 

Here, we would like to generalize the STprior setting to Gabor spaces, which we also term in this context 
shift-and-modulation-invariant (SMI) spaces. We will use these spaces as priors on our input signals, in order to 
recover them from their non-invertible Gabor transform. An SMI subspace W C L 2 is the set of signals that can 
be represented in the form 

fit) = J2 h k ,iM bl T ak w(t), (72) 
fc,zez 
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for some sequence hk,i in t 2 {l?), where w(t) is an arbitrary window in Sq. In other words, W is the closed linear 
span of the Gabor system Q(w, a, b). Our choice of terminology follows from the fact that if f(t) lies in W, then the 
function M b iT ak f(t) is also an element of W for every fixed k, I G Z. Indeed, let f(t) = Yl m nei hm,nM bn T am w(t) 
for some sequence h m , n , then 



M bl T ak f(t) = M bl T ak 

M bn T a mW(t) 

y m,n£Z 

= Y h rn,nM U T ak M bn T am w(t) 

= h?n,ne~ 2niabkn M^ n +i)T a (, m + k) w(t) 
= h m -k,n-ie- 2mah ^ k M bn T am w{t) 

= Y d m,nM bn T am w{t) e W, (73) 

m,n£Z 

where d m , n = h m ^ n ^e- 2 ^ ab ^ k . The same holds for T ak M bl f{t). 

Our setting is thus as follows. We assume that f(t) lies in some SMI space W, generated by Q(w, a, b), which 
we term the prior space, and that we are given the Gabor coefficients c k ,i of f(t), which were computed with the 
analysis window s(t). Our goal is to produce a recovery f{t) using the synthesis window v(t). Clearly, if W does 
not coincide with our synthesis space V, then the reconstruction f(t) cannot equal f(t). The interesting question 
is whether we can obtain the best possible recovery, which is the orthogonal projection / = Pyf, from the Gabor 
coefficients c k j of f(t). As above, we discuss the integer and rational under-sampling cases separately. 

A. Integer under-sampling 

As discussed in Section IIV-CI if the analysis and prior spaces satisfy S © W = L 2 (R), then the recovery 
/ = P v f can be generated by applying the operator H suh = (V*V)~ 1 V*W(S*W)~~ 1 on the Gabor coefficients 
Ck,i prior to synthesis. From Proposition IV. II we know that this direct-sum condition is satisfied if and only if 
& sw (u>,x) 7^ for all u> and x, where <& sw (lo : x) is as in ( |48T > with v(t) replaced by w(t). The operators V*V, 
V*W and S*W correspond to 2D convolutions with the kernels r vv [k, I], r vw [k, 1} and r ww [k, I] respectively, which 
are given by (l49b with the appropriate substitution of s(t), v(t) and w(t). Hence, the operator _ff su b corresponds to 
2D convolution with the filter h su b, whose 2D DTFT is given by 

x ) = ^7 v (74) 

<P SW {ui, X)<P VV [UJ, X) 

where <$> vw (lu, x), <$> sw (uj,x), and <$> vv (uj,x) are the 2D DTFTs of r vw [k, l], r sw [k,l] and r vv [k,l] respectively. 

When the synthesis space V coincides with the prior space W, we have H su \, = (V*V)~ 1 V*W(S*W)~ 1 — 
(S*W)~ 1 , so that the correction filter is the same as in the consistency approach of section \V^A\ In this case, the 
direct-sum condition (namely the invertibility of the operator S*W) guarantees perfect recovery of f(t). To see 
this, note that any / € W can be written as / = Wd for some sequence d^i, so that the Gabor coefficients c& j are 
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given by c = S* f = S*Wd. Therefore, the expansion coefficients can be perfectly recovered using d = (S*W)~ 1 c. 
This property is, of course, independent of the sampling lattice and holds true also in the rational under-sampling 
regime. 

B. Rational under- sampling 

We now extend the subspace-prior approach to the rational under-sampling regime. As before, we assume that the 
input f(t) can be expressed in the form (l72t for some sequence hk,i, where w(t) is a given window in So- As we have 
seen, the best possible recovery, which is the orthogonal projection / = Pyf, can be obtained if the analysis and prior 
spaces satisfy <S ®W = L 2 (M), which in our case is equivalent to r sw [k, I] being invertible with respect to twisted 
convolution. In this case, f = Pvf can be produced by applying the operator i? su b = (V*V)~ 1 V*W(S*W)~~ 1 
on the Gabor transform Ck,i prior to reconstruction. The operators V*V, V*W, and S*W correspond to twisted 
convolution with the kernels r vv [k, I], r v>w [k,l] and r s , w [k,l] respectively. Therefore, H su b corresponds to twisted 
convolution with 

h su b[k, 1} = (r~ x \ r vw \ r~2) [k, I], (75) 

where, r~J-[k,l] and r~^[k,l] are the inverses of r vv [k,l] and r sw [k,l] with respect to \\. Consequently, the $ 
function of the subspace-prior filter is given by 

H mx (cu, x) = 3> sw {lu, x)- 1 ®™^, x)$ vv (lj, x)-\ (76) 

where $ sw (u, x), $ vw (li,x), and & vv (lu,x) are the ^-representations of r sw [k,l], r vw [k,l] and r vv [k,l] respec- 
tively. 

VIII. Twisted Convolution 

In the previous sections, we saw that in order to process the samples c m ,n one needs the inverse of certain 
cross-correlation sequences with respect to \\. In this section we show how to obtain explicitly the inverse of some 
sequence d^.i with respect to twisted convolution with parameter ab. This depends very much on ab. If ab G N, 
then the twisted convolution is a standard convolution, and the Fourier transform can be used to compute the inverse 
of dk y i- If ab = q/p, then one can use the construction derived in [27], which breaks the problem into computing 
inverses of several sequences with respect to standard convolution. We now briefly review this method. For the 
proofs and more detailed explanations, we refer the reader to the original paper. 

Let dk,i be a sequence in ^ 1 (Z 2 ). We create p 2 new sequences out of dk,i, defined as 

{d r ' s )kd = d k j ^ ~ r ~P m i l ~ s ~pn], (77) 

where r, s = 0, 1, . . . ,p — 1. It is easy to see that the sequence d r ' s is supported on the coset (r +p1i) x (s 
and therefore d = YX=0 E!=o ^ r,s - ^ n me case when p = 2, out of a sequence dk,i we obtain four subsequences: 
d°<° which is supported on 2Z x 2Z, d - 1 supported on 2Z x (2Z + 1), ef 1 * supported on (2Z + 1) x 2Z and d 1 ' 1 
supported on (2Z + 1) x (2Z + 1). 
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Next, we associate with the sequence dkj a p x p matrix D whose entries are sequences in t 1 : 

p-i 

D r>s = ^ d m > r - s e- 2 ™ nsq/p , (78) 

m=0 

where r — s should be interpreted as modulo p. This matrix is an element of an algebra M. of p x p matrices with 
multiplication of two matrices D and E given by 

p-i 

[D © E) r<s = D T ,i * E ltS , (79) 

m=0 

where * is a standard convolution. It was shown in [27| that an algebra of such matrices is closed under taking 
inverses, meaning that if D is invertible in Ai then its inverse is also an element of M. and its entries are also 
coming from some sequence in £ 1 (Z 2 ). For example, when p = 2 the above matrix takes the form 

(f,o + d i,o d o,i_ d hi 
d o,i + d i,i d o,o_ d i,o 

where we used the fact that since p = 2, q must be odd, and thus e 27Tlms i/ 2 for m, s = 0, 1 takes the values 1 and 
— 1. Note that summing up the elements of the first column gives us back the sequence d. 

It was shown in [27 1 that the invertibility of the sequence dk.i with respect to \ is equivalent to the invertibility 
of the matrix D in this new matrix algebra, which in turn is equivalent to the invertibility of det(£>) in £ 1 (Z 2 , *). 
Therefore, if D is invertible, its inverse can be computed using Cramer's Rule. That is the (r, s) entry of D^ 1 is 
given by 

(D-% s = (dctp))- 1 * det(D(s,r)), (80) 

where D(s,r) is a p x p matrix obtained from D by substituting the sth row of D with a vector of zeros having 5 
on the rth position, and the rth column with a column of zeros having S on the sth position. Note that det(D) is 
a sequence and its inverse in (l80l > is taken with respect to standard convolution. For example, when p = 2 we get 



D = 




(81) 



Thus, 



D- 1 = (dctD)- 1 * 



d o,o_ d tfi _ d o.i +d i,i 
-d^-d 1 ' 1 d^ + d 1 ' 
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p 2 branches: p 4 branches: 

r, s = 0, . . . ,p — 1 for each r, s, 

u.v — 0, . . . ,p — 1 



Fig. 3: A filter-bank realization of twisted convolution between Ck.i and dk\. 

where det(D) = (d°>° + d 1 ' ) * (d°'° - d 1 ' ) - (d - 1 + d 1 ' 1 ) * (d ' 1 - d 1 ' 1 ). Since the matrix algebra M is closed 
under taking inverses, summing up the elements of the first column of D^ 1 results in some sequence e k ,i which 
is the inverse of dk,i with respect to twisted convolution. Therefore, it is enough to compute only this column 
and sum up its entries to get d _1 . In our example with p = 2, the twisted-convolutional-inverse of d equals 
(det(D))- 1 * (d°'° - d 1 ' - d ' 1 - d 1 ' 1 ). 

We mentioned in the previous sections that it is possible to realize twisted convolution with a rational parameter 
ab using a filter bank. Indeed, using the decomposition d77l ) of the sequences, the twisted convolution of two 
sequences c and d, 

(d^c) m ,„ = J2 dKi^n-K n -i^ 2mah{m - k)l - J2 ^id m _ k ^ ie - 2 ™ b ^ k (82) 
fc.zez k,iei 

can be written as 

p— l p— l 

(dt]c)=^ ^ (c r ' s * d«-r,«-«) e -2iri(t.-.)rg/p for u , w = 0, 1, . . . ,p - 1. (83) 
r,s=0 u,v—0 

Therefore, as shown in Fig. [3] each of the p 2 sequences c r ' s , r, s = 0, 1, . . . ,p — 1, is split into p 2 filters associated 
with the sequences d u,v , u,v = 0, 1, . . . ,p — 1. Then, d\\c is obtained by summing over the resulting p 4 output 
sequences. Figure [3] depicts one of the p 4 branches, which corresponds to the indices r, s ,u and v. 

IX. Example: Integer under-sampling with Gaussian windows 

We now demonstrate the prior-free recovery techniques derived in this paper. To retain simplicity we will focus 
on the integer under-sampling scenario. In this regime, the smallest amount of information loss occurs when ab = 2. 
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Therefore, in our simulations we used a = 1 and b = 2. In this setting there are at most half the number of time- 
frequency coefficients for any given frequency range per time unit, than in any invertible Gabor representation. 
Consequently, algorithms operating in the Gabor domain {e.g., for system identification, speech enhancement, blind 
source separation, etc.) will benefit from a reduction of at least a factor of 2 in computational load. On the other 
hand, we expect the norm of the reconstruction error to be roughly on the order of the signal's norm in the worst-case 
scenario, since half of the information is lost in such a representation. 

For tractability, we will work out the case in which the analysis and synthesis are both performed with a Gaussian 
window: 



1 



ill 



sW =7^ exp r^/ <84) 



1 



ill 



" (() =75^r p r^r <85) 

In this scenario, the cross-correlation sequence r sv [k,l] = (v , M& n T am s) , has an analytic expression: 

1 / (ak) 2 + 4ir 2 <j 2 a 2 ibl) 2 \ ( 2nia 2 s abkl\ 

r sv [k,l} = — cxp^- v ' \ vy \exp\ 2T 2 ■ < 86 ) 

y/2n (ct 2 +<t 2 ) { 2(<t 2 +<t 2 ,) J [ a 2 + cr 2 J 

Similarly, r ss [k,l] and r vv [k,l) can be obtained by replacing a s by a v and vice versa. 

The 2D filter h con of d5TI ). corresponding to the consistency requirement, is the convolutional inverse of r ss [k, I]. 
This sequence can be approximated numerically using the discrete Fourier transform (DFT) of the finite-length 
sequence r sv [k,l], (k,l) 6 [—K,K] x [—L,L], for some (large) K and L. To compute the filter haa of d57b . 
corresponding to the minimax-regret approach, we need to invert r ss [k, I] and r vv [k,l], which can be done in a 
similar manner. Note that both h con and h mx are generally complex sequences. Figure H] depicts the modulus \h coa \ 
and \h mx \ for the case a s — 0.1, a v — 2, and ab = 2. 

To see the effect of these two filters, we now examine the recovery of a chirp signal from its non-invertible 
Gabor representation using both methods. Specifically, let 

/(<) = 2cos(i 2 ). (87) 

The Gaussian-window Gabor transform of f(t) has a closed form expression, given by 

1 f ak(ak + 2Mtt) + 2ib 2 l 2 ir 2 a 2 

c k,l =— F = 5 ex P 



y/-2ia 2 s + 1 I i + 2a 2 

1 f -aktak - 2blir) + 2ib 2 l 2 ir 2 a 2 1 
H / a ex P i . n 9 1 > • (88) 

The signal /(t) and the modulus of its Gabor transform, \ck,i\, are shown in Fig. [5] Although Ck.i seems to constitute 
a good time-frequency representation of f(t), it is certainly not suited to play the role of the synthesis expansion 
coefficients dk,i- This can be seen in Fig. |6ja), where Ck.i have been used without modification as expansion 
coefficients to produce a recovery fit). The signal-to-noise ratio (SNR) of this recovery is 20 log 10 (||/||/||/ — /||) = 
-0.44dB. 

The reconstructions obtained with the consistency and minimax-regret methods are shown in Fig. |6jb) and 
(c). Clearly, they both bear better resemblance to fit). The consistent recovery is the unique signal that can be 
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Fig. 4: The 2D correction filters corresponding to the minimax-regret and consistency methods. 




constructed with the synthesis window v(t), whose Gabor transform coincides with Ck,i- This property makes this 
reconstruction desirable in some sense, although its SNR is only — 1.03dB, worse than the uncompensated recovery. 
To guarantee that the error between our recovery f(t) and the original signal f(t) is small, for every possible f(t) 
that could have generated Ck,i, one has to use the minimax regret approach, as shown Fig.[6|c). This reconstruction 
achieves an SNR of O.ldB, and thus is better than the other two methods in terms of reconstruction error. Figure [7] 
depicts the expansion coefficients dk.i corresponding to the two methods. 
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(c) 

Fig. 6: Reconstructions of f(t) from its Gabor coefficients Ck,i- (a) Without processing Ckj. (b) Consistent recovery, 
namely using dkj = (c * h COD )k,i as expansion coefficients, (c) Minimax-regret recovery, namely using d^.i = 
(c * h mx )f.j as expansion coefficients. 
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Fig. 7: The modulus of the expansion coefficients, \dk,i\, corresponding to the consistent and minimax-regret recovery 
methods. 
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X. Conclusions 

In this paper we explored various techniques for recovering a signal from its non-invertible Gabor transform, 
where the under-sampling factor is rational. Specifically, we studied situations where both the analysis and synthesis 
windows of the transform are given, so that the only freedom is in processing the coefficients in the time-frequency 
domain prior to synthesis. We began with the consistency approach, in which the recovered signal is required 
to possess the same Gabor transform as the original signal. We then analyzed a minimax strategy whereby a 
reconstruction with minimal worst case error is sought. Finally, we developed a recovery method yielding the 
minimal possible error when the original signal is known to lie in some given Gabor space. We showed that all three 
techniques amount to performing a 2D twisted convolution operation on the Gabor coefficients prior to synthesis. 
When the under-sampling factor of the transform is an integer, this process reduces to standard convolution. We 
demonstrated our techniques for Gaussian-window transforms in the context of recovering a chirp signal. 

Appendix A 

The Multiplication Property of the $ representation 

Let Ckj and dk,i be two sequences having # matrix-valued function representations <& c and respectively. 
Then the matrix-valued function 3? associated with the twisted convolution c\d, can be expressed as 

$^ d \u,x) = $ d (u,x)Q> c (u,x). (89) 

Indeed, let again ab — q/p and let r, s = 0, . . . ,p — 1 be fixed, then 



$l c } d) {oj,x)= ^{c\\d)[s-r+pk,l]e 



-2niabrl —2iri(blx-\-akuj) 



EV^ r fj , . , p -2iriab(s-r-\-pk-m)n 
/ L m , n ^s—r -\-pk —m,l — n t> 



— 2itiab(s — rA-pk — 7n)n p — 2itiabrl^—2-Ki{blx-\-akuj) 



p-i 



EY 1 j -2niab{s-r-u)n 

/ . / . ^u-\-pm,n (Jl s — r—u+p(k—m),l—nV & 

ii=0 k,l£Zm,nGZ 
p-1 

^ ^ ^ C-s—u+pmjidu 
u=0 k,l£Z m,n& 



2'Kiab{s—r~-u)n f) — 2-Kiabrl^--2-Ki{blx-\-akuj) 



— 27Tiab(u—r)n —2niabr(l-{-n) — 27ri(b(l-\-n)x-\-a(k-\-m)uj) 
r-\-pk,l& & " 



El \ ^ j — 2-niabrl — 2izi(blx+aku)) \ I \ ^ — 27ziabun — 2-ni(bnx+amu)) 

I / J ^u—r-\-pk,L^ & I I / l-s— u-\-pm,n£ ' 



u—0 \k,l£'L J \m,n€ 

p-1 



S*r,u( w . !B )*S,.( w . a! ) 



u=0 

Hence, ^ d) {tu,x) = <t> d {uj,x)<f? c {uj,x). 

Appendix B 
Proof of Proposition IIII.3I 
Since Q(v, a, b) is a Riesz basis for V, there exist bounds A > and B < oo such that AI p < & vv (li, x) < BI p , 
where & vv (uj,x) is the matrix-valued function associated to the sequence r vv [k, I], defined in d24l) . The system 
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G(w,a,b), with w(t) = Y^,k let hk,lMuT a kV(t), is a Riesz basis if and only if there exist constants C > and 
D < oo such that 

CI P < $ ww {uj,x) <DI p , (90) 

where *b w (ui, x) is a matrix-valued function built from the cross-correlation sequence r ww [k,l] — (w , MbiT a kw) . 
By substituting w(t) = J2k,ie7, h k,iM u T ak v(t) in t ww [k, I] one obtains 

r ww [k, 1} = Yl r ™ ^ -m,z- n]^ m: „e- 27rMfc ( z -") m / ly _ fc , z _ ; e 27rMfc ^- i ) fe 

= ^ (^tl/Ofo.zlVM-ie 2 ™ 6 ^ -0 * 



(/i^r w ^)[fc,Z], (91) 



where r vv [m,n] = (v, Mi> n T am v) and h*[k,l] = h_k,-i- It is easy to check, and we leave it for the reader, that 
& h '(w,x) = $ h (cj,x) H . Therefore, using the relation from Appendix lAl the (r, s)-entry of the matrix <f> ww (u>, x) 
is 

*™(w,x)=(* h (w,x)* m (u,x)*\u,x) H ) , (92) 

where & (u>,x) is a matrix-valued function associated to the sequence hk,i and defined in the Proposition. Hence, 
if G(w, a, b) and Q(v, a, b) are Riesz bases with bounds C > 0, D < oo, and A > 0, B < oo respectively then 



C 



* h (u,x)* h (u,x) H < ^<f> h (co,x)3> m '(LO,x)* h (u;,x) H = ^* ww (tJ,x) < £; 



(93) 
(94) 



Therefore & h (w,x) satisfies (EB with bounds m = C/B and M = D/A. 
On the other hand, if the sequence h^i is such that OTT l is satisfied, then 

<f? ww (uj, x) = $ h (uj, x)$ vv {lu, x)$ h {iu, x) H > A<f> h (tu, x)$ h (cj, x) H > Am (95) 
$ ww (lu, x) = <f> h {uj, x)$ vv (u>, x)$ h (cj, x) H < B$ h (cj, x)$ h {u), x) H < BM, (96) 

and so Q(w,a,b) is a Riesz basis with bounds C = Am and D = BM. It remains to show that Q(w,a,b) and 
G(v, a, b) span the same space. Every element of G(w, a. b) can be uniquely represented by a linear combinations 
of the elements from Q(v,a,b), since the latter is a Riesz basis. It suffices to show that v(t) can be written as a 
linear combination of the elements from Q(w,a,b) (it will be a unique representation since Q(w,a,b) is a Riesz 
basis). Then, since Gabor spaces are closed under translation and modulations, other basis elements from Q(v, a, b) 
will also admit a unique representation in terms of Q(w,a,b). Let g^i be the inverse of hf-i with respect to \\, 
meaning h\g = S. The inverse exists because hf.j satisfies (fJTJ. Let v(t) = J2 m nez 9m,nMi, n T arn w(t). We will 
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now show that v(t) =v{t). Indeed, 



V(t) = ^ 9m,nM bn T am w(t) = ^ E 9rn,nh k jM bn T am M b lT a kv(t) 




= E 3m ' n E hk < ie 2mabml M b(n+l)Ta(m+k)V{t) 



ra.n^L k,l£Z 



= E E g m -k,n-ih k je- 2 ^ m -^M bn T am v(t) 



m,n£Z k,l&Z 



= (h\\g)[m, n]M bn T am v(t) = v(t). 



(97) 



m,n£Z 
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